// Локальные
#include "linal/csrmatrix.h"

namespace linal
{
    //============================================================================================================================
    //============================================================================================================================
    void CSRMatrix::makeIdentity(int row)
    {
        m_di[row] = 1.0;

        for(int j = m_igl[row]; j < m_igl[row + 1]; j++)
            m_ggl[j] = 0;

        for(int j = m_igu[row]; j < m_igu[row + 1]; j++)
            m_ggu[m_ngu[j]] = 0;
    }

    //============================================================================================================================
    //============================================================================================================================
    CSRMatrix CSRMatrix::partFactorLU() const
    {
        CSRMatrix result(m_Dimension, m_igl, m_jgl, m_igu, m_ngu);
        for(int k = 0; k < m_Dimension; k++)
        {
            for(int j = m_igl[k]; j < m_igl[k + 1]; j++)
            {
                result.m_ggl[j] = m_ggl[j];
                result.m_ggu[j] = m_ggu[j];

                for(int i = m_igl[k]; i < j; i++)
                    for(int m = m_igl[m_jgl[j]]; m < m_igl[m_jgl[j] + 1]; m++)
                        if(m_jgl[i] == m_jgl[m])
                        {
                            result.m_ggl[j] -= result.m_ggl[i]*result.m_ggu[m];
                            result.m_ggu[j] -= result.m_ggu[i]*result.m_ggl[m];
                            m = m_igl[m_jgl[j] + 1];
                        }
                        else
                            if(m_jgl[i] < m_jgl[m])
                                m = m_igl[m_jgl[j] + 1];

                result.m_ggu[j] /= result.m_di[m_jgl[j]];
            }

            result.m_di[k] = m_di[k];
            for(int i = m_igl[k]; i < m_igl[k + 1]; i++)
                result.m_di[k] -= result.m_ggl[i]*result.m_ggu[i];
        }
        return result;
    }

    //============================================================================================================================
    //============================================================================================================================
    double& CSRMatrix::operator()(int i, int j)
    {
        if(i == j)
            return m_di[i];
        else if(i > j)
        {
            int index = m_igl[i];
            while(m_jgl[index] != j && index < m_igl[i + 1])
                index++;
            Q_ASSERT_X(index != m_igl[i + 1], "CSRMatrix::operator()", "the element hasn't been found");
            return m_ggl[index];
        }
        else
        {
            int index = m_igl[j];
            while(m_jgl[index] != i && index < m_igl[j + 1])
                index++;
            Q_ASSERT_X(index != m_igl[j + 1], "CSRMatrix::operator()", "the element hasn't been found");
            return m_ggu[index];
        }
    }

    //============================================================================================================================
    //============================================================================================================================
    const double& CSRMatrix::operator()(int i, int j) const
    {
        if(i == j)
            return m_di[i];
        else if(i > j)
        {
            int index = m_igl[i];
            while(m_jgl[index] != j && index < m_igl[i + 1])
                index++;
            Q_ASSERT_X(index != m_igl[i + 1], "CSRMatrix::operator()", "the element hasn't been found");
            return m_ggl[index];
        }
        else
        {
            int index = m_igl[j];
            while(m_jgl[index] != i && index < m_igl[j + 1])
                index++;
            Q_ASSERT_X(index != m_igl[j + 1], "CSRMatrix::operator()", "the element hasn't been found");
            return m_ggu[index];
        }
    }

    //============================================================================================================================
    //============================================================================================================================
    CORE_EXPORT const Vector operator*(const CSRMatrix &matrix, const Vector &vector)
    {
        Q_ASSERT_X(matrix.dimension() == vector.dimension(), "operator*", "matrix and vector have different dimensions");

        Vector result(matrix.dimension());
        for(int i = 0; i < result.dimension(); i++)
            result[i] = vector[i]*matrix.di()[i];

        for(int i = 0; i < result.dimension(); i++)
            for(int j = matrix.ig()[i]; j < matrix.ig()[i + 1]; j++)
            {
                result[i] += vector[matrix.jg()[j]]*matrix.ggl()[j];
                result[matrix.jg()[j]] += vector[i]*matrix.ggu()[j];
            }

        return result;
    }
}
